Variational finite-difference representation of the kinetic energy operator 
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A potential disadvantage of real-space-grid electronic structure methods is the lack of a variational 
principle and the concomitant increase of total energy with grid refinement. We show that the 
origin of this feature is the systematic underestimation of the kinetic energy by the finite difference 
representation of the Laplacian operator. We present an alternative representation that provides 
a rigorous upper bound estimate of the true kinetic energy and we illustrate its properties with 
a harmonic oscillator potential. For a more realistic application, we study the convergence of the 
total energy of bulk silicon using a real-space-grid density-functional code and employing both the 
conventional and the alternative representations of the kinetic energy operator. 



Electronic structure methods based on finite differ- 
ences on a real space grid have gained much support in 
recent years pi |) due to their simplicity and versatil- 
ity. As with plane wave basis sets, their accuracy can be 
improved easily and systematically. In fact, there exists 
a rigorous cutoff for the plane waves that can be repre- 
sented in a given grid, without aliasing, that provides a 
convenient connection between the two schemes. Soft || 
and ultrasoft pseudopotentials, developed in the plane 
wave context, can be applied equally well in grid-based 
methods, resulting in an accurate and efficient evaluation 
of the potential energy. In contrast with plane waves, 
the evaluation of the kinetic energy by finite differences 
is approximate, but it can be significantly improved by 
using high order representations of the Laplacian opera- 
tor (j, 0. However, an important difference between 
finite difference schemes and basis set approaches is the 
lack of a Rayleigh-Ritz variational principle in the former 
case. With finite differences, the accuracy of the calcu- 
lation can also be improved systematically by increasing 
the grid cutoff (i.e. the grid density). But denser grids 
generally result in higher energies. This feature, common 
to all existing real-space-grid approaches to electronic 
structure calculations, has been discussed frequently in 
the literature (see for example the discussion of equation 
(20) in Ref. (|]). It is one of the reasons why it is diffi- 
cult to develop extrapolation methods and convergence 
schemes based on minimizing the total energy. 

In this paper we show that the origin of the "anti vari- 
ational" behavior of the total energy in real-space-grid 
approaches lies in a systematic underestimation of the 
kinetic energy by the finite-difference representation of 
the Laplacian operator, independently of the order used. 
We propose a simple way to construct, for any order, 
an alternative Laplacian representation that leads to ki- 
netic energies that are higher or equal to the true ki- 
netic energy. The paper is organized as follows. We first 
briefly present a way to construct the conventional finite 
difference representation of a Laplacian and show that 
the spectrum of the resulting operator is lower than the 



true one. We then apply a variation of this theme to 
construct a representation which has a spectrum higher 
than the true one. We compare the convergence of the 
one-dimensional harmonic oscillator energy using the two 
representations. Finally we implement the new represen- 
tation in a real-space electronic structure code and ex- 
amine the convergence properties of the calculations of 
bulk Si. 

We will consider a three-dimensional (3D) Cartesian 
coordinate system, in which the Laplacian is the sum of 
three one-dimensional (ID) operators. For a regular grid 
in one dimension, a general finite-difference expression 
for the Laplacian of function ^(cc) at point Xi is 

1 N 

a 3=-N 

where tpi = ip(xi). The constant 1/et 2 takes care of the 
dependence on the grid interval a, so that the coefficients 
Cj are independent of it. Thus, we will use a = 1 for 
simplicity in what follows. 

One way to obtain the Laplacian coefficients Cj is 
through its eigenvalue equation E(k) = k 2 , where E{k) 
is (except for a constant factor) the kinetic energy of a 
single plane wave: 

N 

E[k) = ~ e - ikx V 2 e ikx = -c - c 3 cos (kj), (2) 

where we have used C-j — Cj due to the parity of the 
Laplacian. Notice that E(k) is periodic by construction, 
with period 2n, and that its slope is necessarily zero at 
the grid's Nyquist limit k — ±7r. Since we have only 
N + 1 discrete coefficients to impose E(k) — k 2 in the 
continuous range [— ir, tt], this can be done only approx- 
imately. In an electronic structure calculation, most of 
the spectral weight is concentrated at low values of k. 
Therefore, a sensible prescription is to require the Lapla- 
cian to be accurate around k = 0, by requiring that the 
value of E(k) — k 2 and its N first even derivatives be zero 
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at this point. The resulting values of Cj are given in Table 
I of Ref. |j| However, because of its periodic character, 
and its zero slope at k — tt, the resulting E{k) is a lower 
bound of the true kinetic energy k 2 , as is shown in Fig. [l]. 
With increasing discretization order TV, E{k) approaches 
k 2 , but the convergence is always from below. 

A simple variation of the above theme leads to a rep- 
resentation whose corresponding kinetic energy operator 
is an upper bound to the true kinetic energy ||. Instead 
of fitting N derivatives of E(k) at k = 0, we fit only the 
value and the first N — 1 even derivatives at k = and, 
additionally, the value at k = tt. In Table [i] we present 
the resulting coefficients for orders N = 1 — 6. 
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TABLE I: Laplacian expansion coefficients leading to upper 
bound representation of the kinetic energy operator. 



In Fig. |l| we show the kinetic energy of a plane wave 
as a function of its wavevector k, in the range [0, tt]. The 
middle line is the exact kinetic energy, the lower curve is 
the energy of the conventional Laplacian representation 
of order N = 6, and the upper curve is the result of the 
new Laplacian representation of order 6. It is clearly seen 
that the new representation gives an upper bound to the 
true kinetic energy. 

The simplest example that clearly demonstrates the 
properties of the new Laplacian in action is the ID har- 
monic oscillator. We show the gradual convergence of the 
calculation using the two representations by increasing 
the number of points used to sample the harmonic oscil- 
lator potential \x 2 . In Fig. || we plot the lowest eigen- 
value, whose converged value is 0.5, versus the number 
of points that sample the interval [-5, 5]. We use a sixth 
order representation of the Laplacian for both calcula- 
tions. The convergence is faster with the conventional 
representation, because it samples better the low energy 
part of the spectrum, but the convergence is from below, 
as expected. In contrast, the new representation con- 
verges from above with increasing hamiltonian size, like 
a basis-set expansion. 

Since the 3D Laplacian in Cartesian coordinates is just 
a sum of three independent ID Laplacians, its represen- 
tation in a uniform grid has 67V + 1 nonzero elements 



FIG. 1: Middle (dashed) line: exact kinetic energy of a plane 
wave E(k) = k 2 . Lower (dot-dashed) line: kinetic energy ob- 
tained using the conventional finite difference Laplacian rep- 
resentation of order N = 6. Upper (solid) line: energy of the 
new Laplacian representation with N = 6. 
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FIG. 2: Convergence of the lowest eigenvalue of the harmonic 
oscillator with the number of grid points. The filled circles 
are the results of the upper bound Laplacian representation, 
the open circles are the results of the conventionalLaplacian. 



in a cross orientation, with 3cn at the center and the TV 
elements along the positive and negative sides of the 3 
axes. We have implemented both the conventional and 
the upper bound representation in a real space electronic 
structure code |J. We use the new representation in the 
calculation of the kinetic energy part of the operator, and 
the conventional representation for the solution of the 
Poisson problem. The Poisson part of the problem typi- 
cally converges from above using conventional real-space 
Laplacians, since the contribution of the plane waves en- 
ters with a prefactor of 1/G 2 , where G is the wavevector 
of the plane wave. For the same reason one can argue that 
the Poisson problem needs a lower discretization in its so- 
lution, since the high energy components are damped by 
1/G 2 . 

In Fig. ^ we present the convergence characteristics 
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FIG. 3: Convergence of the total energy calculations for Sili- 
con with respect to the number of grid points in each dimen- 
sion. The curves, from top to bottom correspond to calcu- 
lations using the new kinetic energy formula (filled symbols) 
with orders 4 (squares), 6 (circles), and then the conventional 
formula (open symbols) of order 6 (circles), 4 (squares), and 
2 (diamonds). The lines are only used as a guide to the eye. 



of a total energy calculation for bulk Si, using norm- 
conserving pseudopotentials fl0[| . An 8 atom diamond 
lattice cubic unit cell is used with a 6 x 6 x 6 fc-point grid 

The lattice constant 



in the Monkhorst-Pack scheme 11 



is fixed at 10.264 a.u. and the real-space grid discretiza- 
tion is gradually increased from 14 to 30 points in each 
dimension: this corresponds to an effective plane-wave 
energy cutoff varying from 18 to 84 Ryd. For simplicity, 
we use the same order for the Laplacian in the kinetic and 
Hartree energies. The lowest curve (left pointing trian- 
gles) shows the results of the second order Laplacian that 
is used in HARES 0. The total energy converges from 
below, implying that the major source of error comes 
from the kinetic energy operator. The two curves imme- 
diately above it (up-pointing triangles and rhombuses) 
show the total energy convergence using the conventional 
formulas of order 4 and 6. The energy converges from be- 
low and quickly reaches its correct value. This is partially 
due to a cancellation of errors between the kinetic and 
Hartree energies, but the convergence is not necessarily 
monotonic as is manifest at 16 and 18 grid points per 
dimension. The uppermost curve (circles) shows the re- 
sults using the new Laplacian of order 4 for the kinetic 
energy, and the conventional Laplacian for the Hartree 
energy. The total energy convergence is slower, similar 
to the second order conventional representation, because 
there is no error cancellation, since both the kinetic and 



Hartree energies converge from above. Also, fewer terms 
are used to approximate the low-fc region, while fixing the 
Laplacian value at k = n distorts considerably the kinetic 
energy eigenspectrum of Fig. [l] at low order expansions. 
Finally, the curve immediately below it (squares) shows 
the total energy convergence for the new representation 
of order 6, which has a very satisfactory numerical con- 
verge from above. 

In summary, we derived a new finite difference rep- 
resentation of the Laplacian that is guaranteed to give 
an upper bound value for the kinetic energy operator. 
We demonstrated this basic property of the ID Lapla- 
cian through the harmonic oscillator potential and im- 
plemented and tested the 3D version in a real-space elec- 
tronic structure code. We find that for sufficiently high 
order of the Laplacian, which for the case of bulk Si is 
6, the new form exhibits very satisfactory numerical con- 
vergence of the total energy, while establishing the vari- 
ational character of the real space method. 
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